"""
CBMPy: Constraint Based Modelling in Python (http://pysces.sourceforge.net/cbm)
============
Copyright (C) 2009-2012 Brett G. Olivier, VU University Amsterdam, Amsterdam, The Netherlands

This program is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.

This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU General Public License for more details.

You should have received a copy of the GNU General Public License
along with this program.  If not, see <http://www.gnu.org/licenses/>

Author: Brett G. Olivier
Contact email: bgoli@users.sourceforge.net
Last edit: $Author: bgoli $ ($Id: CBGLPK.py 990 2012-02-15 13:52:11Z bgoli $)

"""

import os, time
import numpy
import glpk # removed for sphinx
import CBWrite, CBTools
from CBConfig import __CBCONFIG__ as __CBCONFIG__
__DEBUG__ = __CBCONFIG__['DEBUG']
__version__ = __CBCONFIG__['VERSION']


GLPK_STATUS = {
    1 : 'LPS_UNDEF',
    2 : 'LPS_FEAS',
    3 : 'LPS_INFEAS',
    4 : 'LPS_NOFEAS',
    5 : 'LPS_OPT',
    6 : 'LPS_UNBND'}


def glpk_analyzeModel(f, lpFname=None, return_lp_obj=False, with_reduced_costs=True, with_sensitivity=False, del_intermediate=False):
    """
    Optimize a model and add the result of the optimization to the model object 
    (e.g. `reaction.value`, `objectiveFunction.value`) and the stoichiometric 
    matrix is automatically generated. This is a common function available in all solver interfaces.
    
     - *f* an instantiated PySCeSCBM model object
     - *lpFname* [default=None] the name of the intermediate LP file. If not specified '_cplxtmp_.tmp'
     - *return_lp_obj* [default=False] off by default when enabled it returns the GLPK LP object
     - *with_reduced_costs* [default=True] calculate and add reduced cost information to mode
     - *with_sensitivity* [default=False] add solution sensitivity information (not yet implemented)
     - *del_intermediate* [default=False] delete the intermediary files after updating model object, useful for server applications
     
    """
    
    CBTools.addStoichToFBAModel(f)
    fid = f.id
    if lpFname == None:
        ##  lpFname = '_cplxtmp_%s_.tmp' % random.randint(0,100)
        lpFname = '_cplxtmp_.tmp'
    f.id = lpFname
    LPF = CBWrite.WriteModelLP(f, work_dir=os.getcwd())
    LPF = os.path.abspath(LPF)
    flp = glpk_getModelFromLP(LPF)
    f.id = fid
    glpk_Solve(flp)
    glpk_setFBAsolutionToModel(f, flp, with_reduced_costs=with_reduced_costs)
    if del_intermediate:
        os.remove(LPF)
    if return_lp_obj:
        return flp
    else:
        del flp

def glpk_setFBAsolutionToModel(fba, lp, with_reduced_costs=False):
    """
    Sets the FBA solution from a CPLEX solution to an FBA object

     - *fba* and fba object
     - *lp* a GLPK LP object
     - *with_reduced_costs* a boolean [default=False]

    """
    sol, objname, objval = glpk_getOptimalSolution(lp)
    print 'Status:', glpk_getSolutionStatus(lp)
    fba.objectives[fba.activeObjIdx].solution, fba.objectives[fba.activeObjIdx].value = sol, objval
    for r in fba.reactions:
        r.value = sol[r.getPid()]
    if with_reduced_costs:
        RC = glpk_getReducedCosts(lp)
        CBTools.setReducedCosts(fba, RC)
    else:
        CBTools.setReducedCosts(fba, {})

def glpk_getReducedCosts(lp):
    """
    Extract ReducedCosts from LP and return as a dictionary

    - GLPK LP

    """

    s_name = glpk_getLPColNames(lp)
    output = {}
    for v in range(0, len(s_name)):
        output.update({s_name[v] : glpk.glp_get_col_dual(lp.lp, v+1)})
    del s_name
    return output

def glpk_getModelFromLP(lptFile, Dir=None):
    """
    Load a LPT (CPLEX format) file and return a GLPK model
    """
    if Dir != None:
        assert os.path.exists(Dir), '\nIncorrect path'
        lptFile = os.path.join(Dir, lptFile)
    lp = glpk.glpk(lptFile)
    return lp

def glpk_getModelFromObj(lpo):
    """
    Return a GLPK object from a FBA model object (via LP file)
    """
    modId = lpo.id
    lpo.id = '_lptmp_'
    LPF = CBWrite.WriteModelLP(lpo, work_dir=os.getcwd())
    lpo.id = modId
    return glpk_getModelFromLP(os.path.split(LPF)[-1], Dir=os.getcwd())

def glpk_Solve(lp):
    """
    Solve the LP and create a status attribute with the solution status
    """
    lp.solve()
    status = glpk_getSolutionStatus(lp)
    print 'Solution status returned as: %s' % status
    lp.solution_status = status
    if status == 'LPS_FEAS' or 'LPS_OPT':
        print '\n++++ Solution ++++'
        print '\t%s = %f\n' % (glpk.glp_get_obj_name(lp.lp), glpk.glp_get_obj_val(lp.lp))
    else:
        print '\n++++ Check solution status ++++\n'

def glpk_getOptimalSolution(lp):
    """
    From a glpk model extract a tuple of solution, ObjFuncName and ObjFuncVal
    """
    FBA_SOLUT = lp.sol.copy()
    OBJF_NAME = glpk.glp_get_obj_name(lp.lp)
    OBJF_VAL = glpk.glp_get_obj_val(lp.lp)
    return FBA_SOLUT, OBJF_NAME, OBJF_VAL

def glpk_getSolutionStatus(lp):
    """
    Returns one of:
    
     - *LPS_OPT*: solution is optimal;
     - *LPS_FEAS*: solution is feasible;
     - *LPS_INFEAS*: solution is infeasible;
     - *LPS_NOFEAS*: problem has no feasible solution;
     - *LPS_UNBND*: problem has unbounded solution;
     - *LPS_UNDEF*: solution is undefined.    

    """
    return GLPK_STATUS[glpk.glp_get_status(lp.lp)]

def glpk_getLPColNames(lp):
    """Return the column names of an LP constraint matrix"""
    ncol_names = []
    ncol = glpk.glp_get_num_cols(lp.lp)
    for c in range(1, ncol+1):
        ncol_names.append(glpk.glp_get_col_name(lp.lp, c))
    return ncol_names

def glpk_getLPRowNames(lp):
    """Return the row names of an LP constraint matrix"""
    nrow_names = []
    nrow = glpk.glp_get_num_rows(lp.lp)
    for r in range(1, nrow+1):
        nrow_names.append(glpk.glp_get_row_name(lp.lp, r))
    return nrow_names

def glpk_getConstraintMatrixInfo(lp):
    """
    Returns the dimensions and labels of the LP constraint matrix.
    """
    ncol = glpk.glp_get_num_cols(lp.lp)
    nrow = glpk.glp_get_num_rows(lp.lp)
    ncol_names = []
    nrow_names = []
    for c in range(1, ncol+1):
        ncol_names.append(glpk.glp_get_col_name(lp.lp, c))
    for r in range(1, nrow+1):
        nrow_names.append(glpk.glp_get_row_name(lp.lp, r))
    return ncol, ncol_names, nrow, nrow_names

def glpk_writeLPtoLPTfile(lp, lptfname, title='', Dir=os.getcwd()):
    """
    Write out a glpk model as an LPT (CPLEX format) file
    """
    glpk.glp_set_prob_name(lp.lp, title)
    glpk.glp_write_lp (lp.lp, None, os.path.join(Dir, lptfname))
    print 'LPT output as %s' % os.path.join(Dir, lptfname)

def glpk_addOptimalSolutionAsConstraint(lp, tol=None, bound='lower', optPercentage=100.0, reset_of=True):
    FBA_SOLUT, OBJF_NAME, OBJF_VAL = glpk_getOptimalSolution(lp)
    ncol, ncol_names, nrow, nrow_names = glpk_getConstraintMatrixInfo(lp)
    nr_idx = glpk.glp_add_rows(lp.lp, 1)
    if __DEBUG__: print 'new row index', nr_idx
    glpk.glp_set_row_name(lp.lp, nr_idx, OBJF_NAME)
    if tol != None:
        #BEGIN# this code is adapted from the COBRA toolbox (add ref), start
        #~ ub =  numpy.ceil(OBJF_VAL/tol)*tol*optPercentage/100.0
        ub =  numpy.ceil(OBJF_VAL/tol)*tol*optPercentage/100.0
        #~ lb =  numpy.floor(OBJF_VAL/tol)*tol*optPercentage/100.0
        lb =  numpy.floor(OBJF_VAL/tol)*tol*optPercentage/100.0   
        #END# this code is adapted from the COBRA toolbox
        if bound == 'upper':
            glpk.glp_set_row_bnds(lp.lp, nr_idx, glpk.GLP_UP, lb, ub)
        else:
            glpk.glp_set_row_bnds(lp.lp, nr_idx, glpk.GLP_LO, lb, ub)
        ##  if ub > lb:
            ##  glpk.glp_set_row_bnds(lp.lp, nr_idx, glpk.GLP_DB, lb, ub)
        ##  else:
            ##  glpk.glp_set_row_bnds(lp.lp, nr_idx, glpk.GLP_LO, lb, ub)
    else:
        ub = OBJF_VAL*(optPercentage/100.0)
        lb = OBJF_VAL*(optPercentage/100.0) 
        if bound == 'upper':
            glpk.glp_set_row_bnds(lp.lp, nr_idx, glpk.GLP_UP, lb, ub)
        elif bound == 'lower':
            glpk.glp_set_row_bnds(lp.lp, nr_idx, glpk.GLP_LO, lb, ub)
        else:
            glpk.glp_set_row_bnds(lp.lp, nr_idx, glpk.GLP_FX, lb, ub)


    print '\nOBJ_FUNC_INTERVAL: %.9f <= %.9f <= %.9f\n' % (lb, OBJF_VAL, ub)

    new_row_range2 = glpk.intArray(ncol+1)
    for i in range(1,ncol+1):
        new_row_range2[i] = i

    # multiple column objective
    new_row2 = glpk.doubleArray(ncol+1)
    for c in range(1, ncol+1):
        oc = glpk.glp_get_obj_coef(lp.lp, c)
        if oc != 0.0:
            new_row2[c] = oc
        else:
            new_row2[c] = 0.0
            
    # single column objective
    ##  new_row2 = glpk.doubleArray(ncol+1)
    ##  for c in range(1, ncol+1):
        ##  if ncol_names[c-1] == OBJF_NAME:
            ##  new_row2[c] = 1.0
        ##  else:
            ##  new_row2[c] = 0.0

    # update row
    glpk.glp_set_mat_row(lp.lp, nr_idx, ncol, new_row_range2, new_row2)

    # think about this either leave the OF as it is or give it zero coefficients
    if reset_of:
        for c in range(1, ncol+1):
            glpk.glp_set_obj_coef(lp.lp, c, 0.0)
            ##  if ncol_names[c-1] == OBJF_NAME:
                ##  glpk.glp_set_obj_coef(lp.lp, c, 1.0)
            ##  else:
                ##  glpk.glp_set_obj_coef(lp.lp, c, 0.0)
    return lp

def glpk_addConstraint(lp, name, value, bound='lower'):
    ncol, ncol_names, nrow, nrow_names = glpk_getConstraintMatrixInfo(lp)
    nr_idx = glpk.glp_add_rows(lp.lp, 1)
    if __DEBUG__: print 'new row index', nr_idx
    glpk.glp_set_row_name(lp.lp, nr_idx, name)
    if bound == 'upper':
        glpk.glp_set_row_bnds(lp.lp, nr_idx, glpk.GLP_UP, value, value)
    elif bound == 'lower':
        glpk.glp_set_row_bnds(lp.lp, nr_idx, glpk.GLP_LO, value, value)

    new_row_range2 = glpk.intArray(ncol+1)
    for i in range(1,ncol+1):
        new_row_range2[i] = i

    # multiple column objective
    new_row2 = glpk.doubleArray(ncol+1)
    for c in range(1, ncol+1):
        oc = glpk.glp_get_obj_coef(lp.lp, c)
        if oc != 0.0:
            new_row2[c] = oc
        else:
            new_row2[c] = 0.0
            
    # update row
    glpk.glp_set_mat_row(lp.lp, nr_idx, ncol, new_row_range2, new_row2)

    return lp

def glpk_runInputScan(fba, exDict, wDir, input_lb=-10.0, input_ub=0.0, writeHformat=False, rationalLPout=False):
    """
    scans all inputs
    """
    debug_dir = os.path.join(wDir, 'lpt')
    ine_dir = os.path.join(wDir, 'ine')
    rat_dir = os.path.join(wDir, 'eslv')
    if not os.path.exists(debug_dir):
        os.mkdir(debug_dir)
    if writeHformat and not os.path.exists(ine_dir):
        os.mkdir(ine_dir)
    if rationalLPout and not os.path.exists(rat_dir):
        os.mkdir(rat_dir)
        
    optimal_growth_rates = {}
    infeasible_inputs = []
    modname = fba.inputfile_id
    fbaid0 = fba.id
    for input in exDict.keys():
        fba.id = '%s(%s)' % (modname, input)
        ilb = fba.getFluxBoundByReactionID(input, 'lower').value
        iub = fba.getFluxBoundByReactionID(input, 'upper').value
        fba.setReactionBound(input, input_lb, 'lower')
        fba.setReactionBound(input, input_ub, 'upper')
        LPfile = CBWrite.WriteModelLP(fba, work_dir=debug_dir)
        fbalp = glpk_getModelFromLP(LPfile, debug_dir)
        ##  fbalp.solve()
        glpk_Solve(fbalp)
        if fbalp.solution_status == 'GLP_OPT':
            opt = glpk_getOptimalSolution(fbalp)[2]
            if __DEBUG__: print 'Optimal growth rate %s(%s) = %s' % (modname, input, opt)
            optimal_growth_rates.update({input : opt})
            if writeHformat:
                CBWrite.WriteModelHFormatFBA(fba, ine_dir, realnum=False)
                CBWrite.WriteModelHFormatFBA(fba, ine_dir, realnum=True)    
            if rationalLPout:
                tmp = CBWrite.WriteModelLP(fba, work_dir=rat_dir)   
        else:
            print 'GLPK returned an error code: %s' % fbalp.solution_status
            infeasible_inputs.append(input)
            os.remove(LPfile)
        fba.setReactionBound(input, ilb, 'lower')
        fba.setReactionBound(input, iub, 'upper')    
    fba.id = fbaid0
    return optimal_growth_rates, infeasible_inputs

def glpk_runInputScanWithFVA(fba, exDict, wDir, input_lb=-10.0, input_ub=0.0, writeHformat=False, rationalLPout=False, NaNonFail=False):
    """
    scans all inputs
    """
    debug_dir = os.path.join(wDir, 'lpt')
    ine_dir = os.path.join(wDir, 'ine')
    rat_dir = os.path.join(wDir, 'eslv')
    if not os.path.exists(debug_dir):
        os.mkdir(debug_dir)
    if writeHformat and not os.path.exists(ine_dir):
        os.mkdir(ine_dir)
    if rationalLPout and not os.path.exists(rat_dir):
        os.mkdir(rat_dir)

    optimal_growth_rates = {}
    infeasible_inputs = []
    modname = fba.inputfile_id
    fbaid0 = fba.id
    fva_out = []
    for input in exDict.keys():
        fba.id = '%s(%s)' % (modname, input)
        ilb = fba.getFluxBoundByReactionID(input, 'lower').value
        iub = fba.getFluxBoundByReactionID(input, 'upper').value
        fba.setReactionBound(input, input_lb, 'lower')
        fba.setReactionBound(input, input_ub, 'upper')
        LPfile = CBWrite.WriteModelLP(fba, work_dir=debug_dir)
        fbalp = glpk_getModelFromLP(LPfile, debug_dir)
        ##  fbalp.solve()
        glpk_Solve(fbalp)
        if fbalp.solution_status == 'GLP_OPT':
            opt = glpk_getOptimalSolution(fbalp)[2]
            if __DEBUG__: print 'Optimal growth rate %s(%s) = %s' % (modname, input, opt)
            optimal_growth_rates.update({input : opt})
            if writeHformat:
                CBWrite.WriteModelHFormatFBA(fba, ine_dir, realnum=False)
                CBWrite.WriteModelHFormatFBA(fba, ine_dir, realnum=True)    
            if rationalLPout:
                tmp = CBWrite.WriteModelLP(fba, work_dir=rat_dir) 
            glpk_addOptimalSolutionAsConstraint(fbalp, tol=None, bound='lower', optPercentage=100.0, reset_of=True)
            OUTPUT_ARRAY, ncol_names = doFluxVariability(fbalp, ignore_fluxes=[], write_debug_lps=False,\
                                        pause_on_skip=True, model_name=fba.id, work_dir=wDir)
            fva_out.append((input,OUTPUT_ARRAY, ncol_names))
            del OUTPUT_ARRAY, ncol_names
        else:
            print 'GLPK returned an error code: %s' % fbalp.solution_status
            infeasible_inputs.append(input)
            os.remove(LPfile)
        fba.setReactionBound(input, ilb, 'lower')
        fba.setReactionBound(input, iub, 'upper')    
    fba.id = fbaid0
    return optimal_growth_rates, infeasible_inputs, fva_out


# old single column method
"""
def glpk_addOptimalSolutionAsConstraintOld(lp, tol=None, optPercentage=100.0):
    FBA_SOLUT, OBJF_NAME, OBJF_VAL = glpk_getOptimalSolution(lp)
    ncol, ncol_names, nrow, nrow_names = glpk_getConstraintMatrixInfo(lp)
    nr_idx = glpk.glp_add_rows(lp.lp, 1)
    if __DEBUG__: print 'new row index', nr_idx
    glpk.glp_set_row_name(lp.lp, nr_idx, OBJF_NAME)
    if tol != None:
        #BEGIN# this code is adapted from the COBRA toolbox (add ref), start
        lb =  numpy.floor(OBJF_VAL/tol)*tol*optPercentage/100.0
        ub =  numpy.ceil(OBJF_VAL/tol)*tol*optPercentage/100.0
        #END# this code is adapted from the COBRA toolbox
        glpk.glp_set_row_bnds(lp.lp, nr_idx, glpk.GLP_LO, lb, ub)
        ##  if ub > lb:
            ##  glpk.glp_set_row_bnds(lp.lp, nr_idx, glpk.GLP_DB, lb, ub)
        ##  else:
            ##  glpk.glp_set_row_bnds(lp.lp, nr_idx, glpk.GLP_LO, lb, ub)
    else:
        lb = ub = OBJF_VAL
        glpk.glp_set_row_bnds(lp.lp, nr_idx, glpk.GLP_LO, lb, ub)

    print '\nOBJ_FUNC_INTERVAL: %.9f <= %.9f <= %.9f\n' % (lb, OBJF_VAL, ub)

    new_row_range2 = glpk.intArray(ncol+1)
    for i in range(1,ncol+1):
        new_row_range2[i] = i

    new_row2 = glpk.doubleArray(ncol+1)
    for c in range(1, ncol+1):
        if ncol_names[c-1] == OBJF_NAME:
            new_row2[c] = 1.0
        else:
            new_row2[c] = 0.0

    ##  print new_row2
    ##  for c in range(1, ncol+1):
        ##  print new_row2[c],
    ##  print ' '

    ##  print new_row_range2
    ##  for i in range(1, ncol+1):
        ##  print new_row_range2[i],
    ##  print ' '
    # update row
    glpk.glp_set_mat_row(lp.lp, nr_idx, ncol, new_row_range2, new_row2)

    # think about this either leave the OF as it is or give it a zero coefficient
    for c in range(1, ncol+1):
        glpk.glp_set_obj_coef(lp.lp, c, 0.0)
        ##  if ncol_names[c-1] == OBJF_NAME:
            ##  glpk.glp_set_obj_coef(lp.lp, c, 1.0)
        ##  else:
            ##  glpk.glp_set_obj_coef(lp.lp, c, 0.0)
    return lp
"""

def glpk_setObjective(c, id, expr=None, sense='maximize', reset=True):
    """
    Set a new objective function
    
     - *c* a GLPK LP object
     - *id* the r_id of the flux to be optimized
     - *expr* a list of (coefficient, flux) pairs
     - *sense* 'maximize'/'minimize'
     - *reset* [default=True] reset all objective function coefficients to zero

    
    """
    sense = sense.lower()
    if sense == 'max': sense = 'maximize'
    if sense == 'min': sense = 'minimize'
    if sense in ['maximise', 'minimise']:
        sense = sense.replace('se','ze')
    assert sense in ['maximize', 'minimize'], "\nsense must be ['maximize', 'minimize'] not %s" % sense    
    glpk.glp_set_obj_name(c.lp, id)
    ncol, ncol_names, nrow, nrow_names = glpk_getConstraintMatrixInfo(c)
    if reset:
        for v in range(ncol):
            glpk.glp_set_obj_coef(c.lp, v+1, 0.0)
        
    if expr != None:
        #~ expr = [(e[1], e[0]) for e in expr]
        if __DEBUG__: print expr
        #~ c.objective.set_linear(expr)
        #~ del expr
        for e in expr:
            print e[0], ncol_names.index(e[1])
            glpk.glp_set_obj_coef(c.lp, ncol_names.index(e[1])+1, e[0])

    if sense == 'minimize':
        glpk.glp_set_obj_dir(c.lp, glpk.GLP_MIN)
        if __DEBUG__: print 'Set minimizing'
    else:
        glpk.glp_set_obj_dir(c.lp, glpk.GLP_MAX)
        if __DEBUG__: print 'Set maximizing'

    if __DEBUG__: 
        print 'writing lp to fva_base_objAdd.lp'
        glpk_writeLPtoLPTfile(c, 'fva_base_objAdd.lp', title=None, Dir=os.getcwd())



            #~ if cntr != 0:
                #~ glpk.glp_set_obj_coef(lp.lp, old_n+1, 0.0)
            #~ glpk.glp_set_obj_coef(lp.lp, n+1, 1.0)
            #~ glpk.glp_set_obj_name(lp.lp, ncol_names[n])
            #~ # minimize
            #~ glpk.glp_set_obj_dir(lp.lp, glpk.GLP_MIN)

def glpk_FluxVariabilityAnalysis(fba, selected_reactions=None, pre_opt=True, tol=None, objF2constr=True, rhs_sense='lower', optPercentage=100.0, work_dir=None, debug=False):
    """
    Perform a flux variability analysis on an fba model:

     - *fba* an FBA model object
     - *selected reactions* [default=None] means use all reactions otherwise use the reactions listed here
     - *pre_opt* [default=True] attempt to presolve the FBA and report its results in the ouput, if this is diabled and *objF2constr* is True then the vid/value of the current active objective is used
     - *tol*  [default=None] do not floor/ceiling the objective function constraint, otherwise round of to *tol*
     - *rhs_sense* [default='lower'] means objC >= objVal the inequality to use for the objective constraint can also be *upper* or *equal*
     - *optPercentage* [default=100.0] means the percentage optimal value to use for the RHS of the objective constraint: optimal_value*(optPercentage/100.0)
     - *work_dir* [default=None] the FVA working directory for temporary files default = cwd+fva
     - *debug* [default=False] if True write out all the intermediate FVA LP's into work_dir
     - *objF2constr* [default=True] add the model objective function as a constraint using rhs_sense etc. If
       this is True with pre_opt=False then the id/value of the active objective is used to form the constraint

    and returns an array with columns::

     Reaction, Reduced Costs, Variability Min, Variability Max, abs(Max-Min), MinStatus, MaxStatus

    and a list containing the row names.
    """
    if work_dir == None:
        work_dir = os.path.join(os.getcwd(), 'fva')
    else:
        assert os.path.exists(work_dir), '\nWhat did you think would happen now!'
        work_dir = os.path.join(work_dir, 'fva')
    if not os.path.exists(work_dir):
        os.mkdir(work_dir)
    debug_dir = work_dir
    ##  if debug:
        ##  debug_dir = os.path.join(work_dir,'DEBUG')
        ##  if not os.path.exists(debug_dir):
            ##  os.mkdir(debug_dir)
    cpx = glpk_getModelFromObj(fba)
    OPTIMAL_PRESOLUTION = None
    pre_sol = pre_oid = pre_oval = None
    REDUCED_COSTS = {}
    if pre_opt:
        glpk_Solve(cpx)
        if glpk_getSolutionStatus(cpx) == 'LPS_OPT':
            print 'Valid FVA presolution'
            OPTIMAL_PRESOLUTION = True
            pre_sol, pre_oid, pre_oval = glpk_getOptimalSolution(cpx)
            fba.objectives[fba.activeObjIdx].solution, fba.objectives[fba.activeObjIdx].value = pre_sol, pre_oval
            REDUCED_COSTS = glpk_getReducedCosts(cpx)
        else:
            print 'Invalid FVA presolution'
            OPTIMAL_PRESOLUTION = False
            pre_sol = {}
            for r in fba.reactions:
                pre_sol.update({r.getPid : 0.0})
                r.reduced_cost = 0.0
            pre_oval = 0.0
            pre_oid = 'None'
            raise RuntimeError, '\nFVA presolve failed to optimize this model and cannot continue!'
    else:
        pre_sol = {}
        for r in fba.reactions:
            pre_sol.update({r.getPid() : 0.0})
        if objF2constr:
            pre_oval = fba.objectives[fba.activeObjIdx].value
            pre_oid = fba.objectives[fba.activeObjIdx].getPid()
        else:
            pre_oval = 0.0
            pre_oid = 'None'        
    for r in fba.reactions:
        r.value = pre_sol[r.id]

    print '\n*****\nStage %s complete\n*****\n' % 1

    # generate new constraint from old objective value (use non-zero coefficients)
    #~ LCS = cpx.objective.get_linear()
    #~ LCN = cpx.variables.get_names()
    #~ new_constraint = []
    #~ for c in range(len(LCS)):
        #~ if LCS[c] != 0.0:
            #~ new_constraint.append((LCS[c], LCN[c]))

    print '\n*****\nStage %s complete\n*****\n' % 2

    # set objective constraint
    if objF2constr:
        glpk_addOptimalSolutionAsConstraint(cpx, tol=tol, bound=rhs_sense, optPercentage=optPercentage, reset_of=True)
    
    #~ if rhs_sense == 'equal':
        #~ raise NotImplementedError, 'rhs_sense=\'equal\' not yet implemented in GLPK.'
        #~ # cplx_setSingleConstraint(cpx, 'ObjCstr', expr=new_constraint, sense='E', rhs=pre_oval)
        #~ glpk_addOptimalSolutionAsConstraint(lp, tol=None, bound=rhs_sense, optPercentage=optPercentage, reset_of=True)
    #~ elif rhs_sense == 'upper':
        #~ if tol != None:
            #~ ub =  numpy.ceil(pre_oval/tol)*tol*optPercentage/100.0
        #~ else:
            #~ ub = pre_oval*(optPercentage/100.0)
        #~ # cplx_setSingleConstraint(cpx, 'ObjCstr', expr=new_constraint, sense='L', rhs=ub)
        #~ glpk_addOptimalSolutionAsConstraint(lp, tol=None, bound=rhs_sense, optPercentage=optPercentage, reset_of=True)
    #~ elif rhs_sense == 'lower':
        #~ if tol != None:
            #~ lb =  numpy.floor(pre_oval/tol)*tol*optPercentage/100.0
        #~ else:
            #~ lb = pre_oval*(optPercentage/100.0)
        #~ # cplx_setSingleConstraint(cpx, 'ObjCstr', expr=new_constraint, sense='G', rhs=lb)
        #~ glpk_addOptimalSolutionAsConstraint(lp, tol=None, bound=rhs_sense, optPercentage=optPercentage, reset_of=True)

    if debug:
        glpk_writeLPtoLPTfile(cpx, 'fva_base.lp', title=None, Dir=work_dir)

    print '\n*****\nStage %s complete\n*****\n' % 3

    NUM_FLX = len(fba.reactions)
    print 'Total number of reactions: %s' % NUM_FLX
    if selected_reactions != None:
        for r in selected_reactions:
            assert r in fba.N.col, "\n%s is not a valid reaction name" % r
        NUM_FLX = len(selected_reactions)
        print 'Number of user selected variables: %s' % NUM_FLX
    try:
        OUTPUT_ARRAY = numpy.zeros((NUM_FLX, 7), numpy.double)
    except AttributeError:
        OUTPUT_ARRAY = numpy.zeros((NUM_FLX, 7))
    OUTPUT_NAMES = []

    print '\n*****\nStage %s complete\n*****\n' % 4
    
    #~ ncol, ncol_names, nrow, nrow_names = glpk_getConstraintMatrixInfo(lp)
    
    for Ridx in range(NUM_FLX):
        if selected_reactions != None:
            R = selected_reactions[Ridx]
        else:
            R = fba.N.col[Ridx]
       
        OUTPUT_NAMES.append(R)
        MIN_STAT = MAX_STAT = 0
        # MIN
        glpk.glp_set_prob_name(cpx.lp, 'FVA min %s' % R)
        ##  cplx_setObjective(cpx, R, expr=None, sense='min', reset=True)
        # TODO: bgoli: see whether this also works with 'minimize'
        glpk_setObjective(cpx, 'min%s' % R, [(1, R)], 'min', reset=True)
        ##  cplx_setBounds(c, id, min=None, max=None) # think about this
        glpk_Solve(cpx)
        if glpk_getSolutionStatus(cpx) == 'LPS_OPT':
            min_sol, min_oid, min_oval = glpk_getOptimalSolution(cpx)
            MIN_STAT = 1
        else:
            ##  min_sol, min_oid, min_oval = 0.0
            min_sol = None
            min_oid = None
            min_oval = numpy.NaN
        if debug:
            glpk_writeLPtoLPTfile(cpx, '%smin' % R, title='min%s=%s' % (R,min_oval), Dir=work_dir)
        
        print '\n*****\nStage %s complete\n*****\n' % 5

        # MAX
        glpk.glp_set_prob_name(cpx.lp, 'FVA max %s' % R)
        ##  cplx_setObjective(cpx, R, expr=None, sense='max', reset=True)
        glpk_setObjective(cpx, 'max%s' % R, expr=None, sense='max', reset=False)
        ##  cplx_setBounds(c, id, min=None, max=None) # think about this
        glpk_Solve(cpx)
        if glpk_getSolutionStatus(cpx) == 'LPS_OPT':
            max_sol, max_oid, max_oval = glpk_getOptimalSolution(cpx)
            MAX_STAT = 1
        else:
            ##  max_sol, max_oid, max_oval = 0.0
            max_sol = None
            max_oid = None
            max_oval = numpy.NaN
        if debug:
            glpk_writeLPtoLPTfile(cpx, '%smax' % R, title='max%s=%s' % (R,max_oval), Dir=work_dir)

        OUTPUT_ARRAY[Ridx,0] = pre_sol[R]
        if REDUCED_COSTS.has_key(R):
            OUTPUT_ARRAY[Ridx,1] = REDUCED_COSTS[R]
        OUTPUT_ARRAY[Ridx,2] = min_oval
        OUTPUT_ARRAY[Ridx,3] = max_oval
        OUTPUT_ARRAY[Ridx,4] = abs(max_oval - min_oval)
        OUTPUT_ARRAY[Ridx,5] = MIN_STAT
        OUTPUT_ARRAY[Ridx,6] = MAX_STAT
        if MIN_STAT == 0.0 or MIN_STAT == 0.0:
            print 'Solver fail for reaction'
            print R, OUTPUT_ARRAY[Ridx]
            time.sleep(5)
            
    print '\n*****\nStage %s complete\n*****\n' % 6
    
    print 'Output array has columns:'
    print 'Reaction, Reduced Costs, Variability Min, Variability Max, abs(Max-Min), MinStatus, MaxStatus'
    return OUTPUT_ARRAY, OUTPUT_NAMES
    
    
def glpk_doFluxVariability(lp, ignore_fluxes=[], write_debug_lps=False, pause_on_skip=True, NaNonFail=True, model_name='Model', work_dir=None):
    raise DeprecationWarning, '\nWarning: this method is being replaced with the \"glpk_FluxVariabilityAnalysis(...)\"'
    
    if work_dir == None:
        work_dir = os.path.join(os.getcwd(), 'fva')
    else:
        assert os.path.exists(work_dir), '\nWhat did you think would happen now!'
        work_dir = os.path.join(work_dir, 'fva') 
    if not os.path.exists(work_dir):
        os.mkdir(work_dir)
    ncol, ncol_names, nrow, nrow_names = glpk_getConstraintMatrixInfo(lp)
    try:
        OUTPUT_ARRAY = numpy.zeros((ncol, 3), numpy.double)
    except AttributeError:
        OUTPUT_ARRAY = numpy.zeros((ncol, 3))
    OUTPUT_ARRAY[:,:] = numpy.NaN
        
    FBA_SOLUT, OBJF_NAME, OBJF_VAL = glpk_getOptimalSolution(lp)
    cntr = 0
    old_n = None
    LPSOLVECNT = 0
    for n in range(ncol):
        if ncol_names[n] not in ignore_fluxes:
            if cntr != 0:
                glpk.glp_set_obj_coef(lp.lp, old_n+1, 0.0)
            glpk.glp_set_obj_coef(lp.lp, n+1, 1.0)
            glpk.glp_set_obj_name(lp.lp, ncol_names[n])
            # minimize
            glpk.glp_set_obj_dir(lp.lp, glpk.GLP_MIN)
            NoFeasMinSol = False
            # lp.solve()
            glpk_Solve(lp)
            print ' ** Minimizing %s: LP number %s ** ' % (ncol_names[n], LPSOLVECNT)
            LPSOLVECNT += 1
            if lp.solution_status == 'GLP_OPT':
                MIN = glpk.glp_get_obj_val(lp.lp)
            else:
                NoFeasMinSol = True
                ##  MIN = numpy.NaN
                if NaNonFail:
                    MIN = numpy.NaN
                else:
                    MIN = 0.0
            glpk.glp_set_prob_name(lp.lp, 'Minimize %s (feasible=%s)' % (ncol_names[n], NoFeasMinSol))
            if write_debug_lps:
                glpk.glp_write_lp (lp.lp, None, os.path.join(work_dir, '%s(%s).min.lp' % (model_name, ncol_names[n])))

            # maximize
            glpk.glp_set_obj_dir(lp.lp, glpk.GLP_MAX)
            NoFeasMaxSol = False
            # lp.solve()
            glpk_Solve(lp)
            print ' ** Maximizing %s: LP number %s ** ' % (ncol_names[n], LPSOLVECNT)
            LPSOLVECNT += 1
            if lp.solution_status == 'GLP_OPT':
                MAX = glpk.glp_get_obj_val(lp.lp)
            else:
                NoFeasMaxSol = True
                if NaNonFail:
                    MAX = numpy.NaN
                else:
                    MAX = 0.0
            glpk.glp_set_prob_name(lp.lp, 'Maximize %s (feasible=%s)' % (ncol_names[n], NoFeasMaxSol))
            if write_debug_lps:
                glpk.glp_write_lp (lp.lp, None, os.path.join(work_dir, '%s(%s).max.lp' % (model_name, ncol_names[n])))

            OUTPUT_ARRAY[n,0] = MIN
            OUTPUT_ARRAY[n,1] = FBA_SOLUT[ncol_names[n]]
            OUTPUT_ARRAY[n,2] = MAX
            cntr += 1
            old_n = n
            if NoFeasMinSol or NoFeasMaxSol:
                print 'Infeasible %s (min=%s, max=%s)' % (n, NoFeasMinSol, NoFeasMaxSol)
        else:
            print '\n*****\nSkipping flux: %s\n*****\n' % ncol_names[n]
            OUTPUT_ARRAY[n, 1] = FBA_SOLUT[ncol_names[n]]
            if pause_on_skip: time.sleep(0.2)

    return OUTPUT_ARRAY, ncol_names

def glpk_doFluxVariability_alt(lp, ignore_fluxes=[], write_debug_lps=False, model_name='Model', work_dir=None, pause_on_skip=True, pause_on_fail=False):
    raise DeprecationWarning, '\nWarning: this method is being replaced with the \"glpk_FluxVariabilityAnalysis(...)\"'
    
    print '\n\nNote: this method does not use a reference FBA and outputs an array as [MIN, MAX, abs(MAX-MIN)] where MIN or MAX are infeasible the other is used.\n'
    if work_dir == None:
        work_dir = os.path.join(os.getcwd(), 'fva')
    else:
        assert os.path.exists(work_dir), '\nWhat did you think would happen now!'
        work_dir = os.path.join(work_dir, 'fva') 
    if not os.path.exists(work_dir):
        os.mkdir(work_dir)
    ncol, ncol_names, nrow, nrow_names = glpk_getConstraintMatrixInfo(lp)
    try:
        OUTPUT_ARRAY = numpy.zeros((ncol, 3), numpy.double)
    except AttributeError:
        OUTPUT_ARRAY = numpy.zeros((ncol, 3))
    OUTPUT_ARRAY[:,:] = numpy.NaN

    # remove objective function
    for c in range(1, ncol+1):
        glpk.glp_set_obj_coef(lp.lp, c, 0.0)

    cntr = 0
    old_n = None
    LPSOLVECNT = 0
    for n in range(ncol):
        if ncol_names[n] not in ignore_fluxes:
            if cntr != 0:
                glpk.glp_set_obj_coef(lp.lp, old_n+1, 0.0)
            glpk.glp_set_obj_coef(lp.lp, n+1, 1.0)
            glpk.glp_set_obj_name(lp.lp, ncol_names[n])
            # minimize
            glpk.glp_set_obj_dir(lp.lp, glpk.GLP_MIN)
            NoFeasMinSol = False
            # lp.solve()
            glpk_Solve(lp)
            print ' ** Minimizing %s: LP number %s ** ' % (ncol_names[n], LPSOLVECNT)
            LPSOLVECNT += 1
            if lp.solution_status == 'GLP_OPT':
                MIN = glpk.glp_get_obj_val(lp.lp)
            else:
                NoFeasMinSol = True
                # MIN = numpy.NaN # LP state of mind
                MIN = 0.0
            glpk.glp_set_prob_name(lp.lp, 'Minimize %s (feasible=%s)' % (ncol_names[n], NoFeasMinSol))
            if write_debug_lps:
                glpk.glp_write_lp (lp.lp, None, os.path.join(work_dir, '%s(%s).min.lp' % (model_name, ncol_names[n])))

            # maximize
            glpk.glp_set_obj_dir(lp.lp, glpk.GLP_MAX)
            NoFeasMaxSol = False
            # lp.solve()
            glpk_Solve(lp)
            print ' ** Maximizing %s: LP number %s ** ' % (ncol_names[n], LPSOLVECNT)
            LPSOLVECNT += 1
            if lp.solution_status == 'GLP_OPT':
                MAX = glpk.glp_get_obj_val(lp.lp)
            else:
                NoFeasMaxSol = True
                # MAX = numpy.NaN
                MAX = 0.0
            glpk.glp_set_prob_name(lp.lp, 'Maximize %s (feasible=%s)' % (ncol_names[n], NoFeasMaxSol))
            if write_debug_lps:
                glpk.glp_write_lp (lp.lp, None, os.path.join(work_dir, '%s(%s).max.lp' % (model_name, ncol_names[n])))

            OUTPUT_ARRAY[n,0] = MIN
            OUTPUT_ARRAY[n,1] = MAX
            if MIN != numpy.NaN and MAX != numpy.NaN:
                OUTPUT_ARRAY[n,2] = abs(MAX-MIN)
            elif MIN == numpy.NaN:
                OUTPUT_ARRAY[n,2] = abs(MAX)
            elif MAX == numpy.NaN:
                OUTPUT_ARRAY[n,2] = abs(MIN)
            cntr += 1
            old_n = n
            if NoFeasMinSol or NoFeasMaxSol:
                print 'Infeasible %s (min=%s, max=%s)' % (n, NoFeasMinSol, NoFeasMaxSol)
                if pause_on_fail:
                    raw_input(n)
        else:
            print '\n*****\nSkipping flux: %s\n*****\n' % ncol_names[n]
            ##  OUTPUT_ARRAY[n, 1] = FBA_SOLUT[ncol_names[n]]
            if pause_on_skip: time.sleep(0.2)

    return OUTPUT_ARRAY, ncol_names


def glpk_cleanFluxVariabilityOutput(flux_variability, log10=False, round=None):
    raise NotImplementedError, '\n Needs to be updated to new solver interface'
    m_arch = numpy.MachAr()
    for r in range(flux_variability.shape[0]):
        for c in range(flux_variability.shape[1]):
            if abs(flux_variability[r,c]) < m_arch.eps/100.0:
                flux_variability[r,c] = 0.0
            if flux_variability[r,c] == numpy.NaN or flux_variability[r,c] == 0.0:
                pass
            elif log10:
                if flux_variability[r,c] > 0.0:
                    flux_variability[r,c] = numpy.log10(flux_variability[r,c])
                elif flux_variability[r,c] < 0.0:
                    flux_variability[r,c] = -numpy.log10(abs(flux_variability[r,c]))
            if round != None:
                flux_variability[r,c] = numpy.round(flux_variability[r,c], round)
            ##  elif abs(flux_variability[r,c]) < m_arch.eps:
                ##  flux_variability[r,c] = 0.0
            ##  else:
                ##  flux_variability[r,c] = numpy.log10(abs(flux_variability[r,c]))
            ##  if abs(flux_variability[r,c]) < 1.0e-9:
            ##  if abs(flux_variability[r,c]) < m_arch.eps/10.0:
                ##  flux_variability[r,c] = 0.0
            ##  if flux_variability[r,c] > 0.0:
                ##  flux_variability[r,c] = numpy.log10(flux_variability[r,c])
            ##  elif flux_variability[r,c] < 0.0:
                ##  flux_variability[r,c] = -numpy.log10(abs(flux_variability[r,c]))
    return flux_variability

def glpk_plotFluxVariability(flux_variability, flux_variability_row, ySlice=None, minHeight=None, maxHeight=None, fname=None):
    raise RuntimeError, '\nObsolete see CBplot module'
    assert _HAVE_MATPLOTLIB_, "\nPlotting requires Matplotlib"
    pylab.ion()
    pylab.hold(False)
    l_cntr = 0
    c_width = 0.8
    g_bars = []
    g_bars_lcorner =[]
    vResults = {}
    PLOTLOG = False
    outputNames = []
    ##  for r in range(flux_variability.shape[0]):
        ##  if flux_variability[r,1] == 0.0:
            ##  print 'No variability: %s' % flux_variability_row[r]
        ##  elif flux_variability[r,0] == flux_variability[r,1] and flux_variability[r,1] == flux_variability[r,2]:
            ##  print 'No variability: %s' % flux_variability_row[r]
        ##  else:
            ##  outputNames.append(flux_variability_row[r])
            ##  print flux_variability_row[r],
            ##  print flux_variability[r]
            ##  b_height = abs(flux_variability[r,2]-flux_variability[r,0])
            ##  g_bars.append(matplotlib.pyplot.bar(left=l_cntr, height=b_height,\
               ##  width=c_width, bottom=flux_variability[r,0], log=PLOTLOG))
            ##  g_bars_lcorner.append(l_cntr)
            ##  l_cntr += c_width
            ##  vResults.update({flux_variability_row[r] : flux_variability[r].copy()})
    Ymagic = []
    for r in range(flux_variability.shape[0]):
        HASMIN = False
        HASMAX = False
        fv_min = flux_variability[r,0]
        fv_fba = flux_variability[r,1]
        fv_max = flux_variability[r,2]
        if fv_fba != numpy.NaN:
            if fv_min != numpy.NaN:
                if fv_min < fv_fba:
                    HASMIN = True
            if fv_max != numpy.NaN:
                if fv_max > fv_fba:
                    HASMAX = True
        b_height = 0.0
        b_height1 = 0.0
        b_height2 = 0.0
        if HASMAX:
            b_height1 = fv_max-fv_fba
        if HASMIN:
            b_height2 = fv_fba-fv_min
        b_height = abs(b_height1)+abs(b_height2)

        HCheckMin = False
        HCheckMax = False
        if minHeight == None:
            HCheckMin = True
        elif minHeight != None and b_height >= minHeight:
            HCheckMin = True
        if maxHeight == None:
            HCheckMax = True
        elif maxHeight != None and b_height <= maxHeight:
            HCheckMax = True
        if b_height > 0.0 and HCheckMin and HCheckMax:
            outputNames.append(flux_variability_row[r])
            if HASMIN:
                bottom = fv_min
            else:
                bottom = fv_fba
            Ymagic.append(bottom)
            Ymagic.append(bottom+b_height)
            ##  print 'Bar = (%s,%s)' % (bottom, bottom+b_height)
            g_bars.append(matplotlib.pyplot.bar(left=l_cntr, height=b_height,\
               width=c_width, bottom=bottom, log=PLOTLOG))
            pylab.hold(True)
            g_bars_lcorner.append(l_cntr)
            l_cntr += c_width
            vResults.update({flux_variability_row[r] : flux_variability[r].copy()})
        ##  if flux_variability_row[r] == 'R_GAPD':
            ##  print vResults['R_GAPD']
            ##  print fv_min, fv_fba, fv_max
            ##  print b_height
            ##  raw_input()
    if __DEBUG__: print 'len flux_variability_row', len(flux_variability_row)
    if __DEBUG__: print 'len g_bars', len(g_bars)
    ##  print 'flux_variability.shape', flux_variability.shape
    outputNames = [l.replace('_LPAREN_e_RPAREN_','_e').replace('R_','') for l in outputNames]
    matplotlib.pyplot.xticks(numpy.array(g_bars_lcorner)+(c_width/2.0), outputNames,\
        rotation='vertical', size='xx-small')
    matplotlib.pyplot.title('%s has %i varying fluxes' % (fname, len(g_bars)))
    matplotlib.pyplot.ylabel('Variability')
    yhi = max(Ymagic) + 0.01*max(Ymagic)
    ylow = min(Ymagic) - abs(0.01*min(Ymagic))
    if ySlice != None:
        yhi = abs(ySlice)
        ylow = -abs(ySlice)
    matplotlib.pyplot.ylim(ylow, yhi)
    if __DEBUG__: print 'Plotting y %s --> %s' % (ylow, yhi)
    if fname != None:
        matplotlib.pyplot.savefig(os.path.join(cDir, fname))

